In [1]:
from IPython.core.display import HTML
css_file = '../ipython_notebook_styles/ngcmstyle.css'
HTML(open(css_file, "r").read())
Out[1]:
Everything that will be done is classical; done in terms of Newtonian physics. Particles will be spherical balls, forces given by harmonic oscillator springs. This cartoon model is useful for quite complex situations.
The two lectures will be on
Molecular mechanics is a method by which molecular systems are modelled using classical mechanics.
We want to make models to predict the behaviour of molecules (changes in their structure, and how they interact with other molecules).
Can't do this in the lab thanks to cost, fiddling with underlying "physical" parameters (masses etc), etc.
Underlying model is that the potential energy of a molecule can be described by the force field using the principles of classical mechanics.
For a typical force field, the potential energy of the molecule is determined in terms of the
Typically there's a massive table giving all this information
Essential : Cartesian coordinates, and prescription for the force field.
These interactions are not just within a single molecule; inter-molecular forces are also important, as well as the standard intra-molecular forces.
To formally relate potential energy to the positions of the atmos in our molecule of interest we use
$$ \begin{equation} \vec{F} = m \times \vec{a} \end{equation} $$Using that the acceleration $\vec{a}$ is related to the position $\vec{r}$ as
$$ \begin{align} \vec{a} &= \frac{d \vec{v}}{d t}, \\ \vec{v} &= \frac{d \vec{r}}{d t} \end{align} $$and expressing the force as the (negative) gradient of the potential energy,
$$ \begin{equation} \vec{F} = - \frac{d V}{d \vec{r}}, \end{equation} $$we can combine these results to get
$$ \begin{equation} -\frac{d V}{d \vec{r}_i} = m_i \frac{d^2 \vec{r}_i}{d t^2}. \end{equation} $$Here $i$ is thought of as an index on the atoms making up the molecule.
For simplicity divide the foce into two parts:
$$ \begin{equation} V_{\text{pot}} = \sum_{\text{bonded}} + \sum_{\text{non-bonded}}. \end{equation} $$The bonded terms include bonds, angles, dihedrals (torsions).
The non-bonded terms include permanent electrostatics, disperion and repulsion.
Usually use a harmonic potential (decent approximation):
$$ \begin{align} E &= \sum_{\text{bonds}} k_b (r - r_0)^2 \\ E &= \sum_{\text{angles}} k_{\theta} ( \theta - \theta_0 )^2. \end{align} $$A Morse function (where the bond "breaks" after a certain length, leading to a constant potential above a certain length) would be more accurate for bond length.
The energetic penalty associated with deviating from the equilibrium has the same form for bond lengths and bond angles as shown above.
Changing the $k_{b, \theta}$ parameters will broaden or steepen the slope of the potential well
Cannot use a simple harmonic term for dihedral angles. This is because there's multiple configurations with the same minima, and two parameters to consider.
Simplest model is to use
$$ \begin{equation} E = \sum_{\text{dihedrals}} A \left[ 1 + \cos (n \phi - \phi_0) \right] \end{equation} $$The $A$ parameter controls the amplitude of the potential, whilst the $n$ parameter controls its periodicity, and $\phi_0$ shifts the curve along the rotation angle axis.
Non-bonded energy represents the pair-wise sum of the energies of all possible interacting (non-bonded) atoms $i$ and $j$. Typically use the Lennard-Jones potential:
$$ \begin{equation} V(r) = 4 \epsilon \left[ \left( \frac{\sigma}{r_{ij}} \right)^{12} - \left( \frac{\sigma}{r_{ij}} \right)^{6} \right]. \end{equation} $$Dispersion is always attractive, occurs at short range and rapidly falls off with distance.
Repulsion occurs when the distance between interating atmos becomes slightly less than the sum of their contact radii.
The $\epsilon$ term defines the "stickiness" of the potential. What is the optimum interaction between the two particles? Do they "like" to be near each other?
The $\sigma$ term defines the size of the particles. Hard spheres cannot inter-penetrate, so the point at which the energy begins to rise steeply is the limit of the particle.
There are (large) tables of parameters giving the values of $\sigma$ and $\epsilon$ for given atoms. When the atmos mix you typically use the Lorenz-Berthelot mixing rules:
$$ \begin{align} \sigma_{AB} & = \tfrac{1}{2} \left( \sigma_{AA} + \sigma_{BB} \right), \\ \epsilon_{AB} & = \sqrt{\epsilon_{AA} \epsilon_{AB}}. \end{align} $$Note: you cannot take one part of a force field (eg the bond length) and mix with another part of a different force field (eg the non-bonded terms). They're self-consistent and come as a whole.
Even when a molecule does not have a formal charge, because the charge distributions are not uniform the atoms have different electronegativities. This leads to atoms have partial charges. These charges are additive pairwise interactions again, but could be attractive or repulsive.
Two representations:
The former is easier to code and explain, and more efficient to code. In this case we have
$$ \begin{equation} V = \sum_{\text{pairs } (i, j)} \frac{q_i q_j}{4 \pi \epsilon_0 r_{ij}}. \end{equation} $$This is simply Coulomb's law.
This is often the most costly term in the system to solve due to the long range interactions; often the least interesting as well, but crucial.
Various ways of dealing with the long-range terms, including
Need
Then
gives as output
This then needs to be analysed by the user.
The DEs from Newton's second law are solved.
Trajectory that specifies hwo the positions and velocities of the atoms vary with time.
Equations solved by taking a small step in time, predicting new atom positions and velocities at the end of the step.
At the new positions, re-calculate the atomic forces and take another timestep.
Coupled nature of the problem, so we solve them numerically.
Assume that the functions $\vec{r}, \vec{v}, \vec{a}$ can be approximated using a Taylor series expansion
$$ \begin{align} r(t + \delta t) & = r(t) + v(t) \, \delta t + \tfrac{1}{2} a(t) \, \delta t^2, \\ v(t + \delta t) & = v(t) + a(t) \, \delta t + \tfrac{1}{2} b(t) \, \delta t^2, \\ a(t + \delta t) & = a(t) + b(t) \, \delta t + \dots \end{align} $$The Verlet and velocity Verlet algorithms are the most popular.
Uses the positions and accelerations at time $t$ and positions from time $t - \delta t$ to calculate new positions at time $t + \delta t$.
This does not use explicit velocities.
Algorithm uses the Taylor series expansion for $r$ at $t \pm \delta t$ which, when summed, give
$$ \begin{equation} r(t + \delta t) + r(t - \delta t) = 2 r(t) + a(t) \, \delta t^2. \end{equation} $$Advantages
Disadvantages